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Vortex arrays in type-II superconductors admit the translational symmetry of an 
infinite system. There are cases, however, like ultra-cold trapped Fermi gases and the 
crust of neutron stars, where finite-size effects make it quite more complex to account 
for the geometrical arrangement of vortices. Here, we self-consistently generate these 
arrays of vortices at zero and finite temperature through a microscopic description of 
the non-homogeneous superfluid based on a differential equation for the local order 
parameter, obtained by coarse graining the Bogoliubov-de Gennes (BdG) equations. 
In this way, the strength of the inter-particle interaction is varied along the BCS-BEC 
crossover, from largely overlapping Cooper pairs in the BCS limit to dilute composite 
bosons in the BEC limit. Detailed comparison with two landmark experiments on 
ultra-cold Eermi gases, aimed at revealing the presence of the superfluid phase, brings 
out several features that makes them relevant for other systems in nature as well. 


The Meissner effect provides evidence of the supercon¬ 
ducting phase below the critical temperature Tc [U |2] . In 
particular, in type-II superconductors the flux of an ap¬ 
plied magnetic held H is partially or totally expelled from 
the material, depending on its strength. Two critical val¬ 
ues of H are distinguished, such that for 0 < i/ < iJci the 
magnetic flux is totally expelled while for 11^ < H < Hc 2 
the system allows for the penetration of an ordered array 
of vortices. Both and depend on the tempera¬ 
ture T, making Tc depend on H. 

These considerations apply to large systems contain¬ 
ing an (essentially) inhnite number of particles, such 
that also the number of vortices can be considered in¬ 
finite for all practical purposes. There exist, however, 
systems with a finite number of particles N where su- 
perfiuid behavior occurs. They include nuclei [3] with 
N ~ 10^ and ultra-cold trapped Fermi gases m typically 
with N ~ 10"^ —10®. In these systems, orbital effects, 
which in superconducting materials are associated with 
an applied magnetic field, result by rotating the sample 
as a whole [5] or through synthetic gauge fields [SHE]. 

In particular, ultra-cold Fermi gases have attracted 
special interest recently [9l [10], because their inter¬ 
particle interaction can be varied almost at will via the 
use of Fano-Feshbach resonances [IIHIlj- These cir¬ 
cumstances make the Meissner-like effect in these sys¬ 
tems depend also on the coupling parameter 
(see Methods). In addition, surface effects may become 
prominent in systems with finite N, since a finite fraction 
of the sample can become normal in the “outer” region 
due to rotation while the “inner” region remains super- 
fluid. The Meissner-like effect in these bounded systems 
can be connected with related effects occurring in nuclei 
where the moment of inertia gets quenched by superflu¬ 
idity [3], and in the inner crust of neutron stars where 
array of vortices form close to the surface [ISlilS]. 

Two experiments were performed by setting an ultra¬ 


cold trapped Fermi gas into rotation, aiming at revealing 
the presence of the superffuid phase. In the first experi¬ 
ment arrays of vortices were generated at low tem¬ 
perature throughout the BCS-BEC crossover, providing 
the first direct evidence for the occurrence of the super¬ 
fluid phase in these systems. In the second experiment 
[18] , the superfluid behavior was revealed by the quench¬ 
ing of the moment of inertia, which was measured at uni- 
tarity (where {kpap)~^ = 0) for increasing temperature 
until its classical value of the normal state was reached. 
In both experiments, the two lowest hyperfine atomic lev¬ 
els were equally populated for a gas of (®Li) Fermi atoms 
contained in an ellipsoidal trap of radial Vr and axial 
frequencies (see Methods), and the system was assumed 
to reach a (quasi) equilibrium situation at any given ro¬ 
tation frequency iz. These experiments complement each 
other, in that they involve determining the total angular 
momentum and controlling the appearance of vortices. 

Despite the role played by these experiments, no theo¬ 
retical account exists for the structure of complex arrays 
of quantum vortices in a rotating superfluid with a fi¬ 
nite number of particles in realistic traps as functions of 
inter-particle coupling and temperature, nor for the re¬ 
lated behavior of the moment of inertia. This is because 
accounting theoretically for these complex arrays of vor¬ 
tices across the BCS-BEC crossover becomes computa¬ 
tionally prohibitive when relying on a standard tool like 
the BdG equations, which extend the BCS mean field to 
nonuniform situations [T9|. Here, we tackle this difficult 
computational problem by resorting to a novel differen¬ 
tial equation for the local gap parameter A(r) at posi¬ 
tion r, through which it is possible to find solutions in 
a rotating trap also for large numbers of vortices. This 
equation was introduced in ref.[3D] (see Methods) and 
rests on the assumption that the phase (p{r) of the com¬ 
plex function A(r) = | A(r)| varies more rapidly 

than the magnitude |A(r)|. Accordingly, this equation 
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was referred to as a Local Phase Density Approximation 
(LPDA) to the BdG equations. In the following, we shall 
determine the behavior of the gap parameter, density, 
and current for a rotating trapped Fermi gas using the 
LPDA equation, while varying the thermodynamic vari¬ 
ables {kpaF)~^ and T, as well as the angular frequency 
D = 2Trh'. 

Our theoretical analysis highlights several features 
which are of interest also to other branches of physics, 
such as: (i) The deviations from Feynman’s theorem for 
the vortex spacing in an array; (ii) The bending of a vor¬ 
tex filament at the surface of the cloud; (hi) The shape of 
the central vortex surrounded by a finite number of vor¬ 
tices; (iv) The emergence of a bi-modal distribution for 
the density profile in the superfluid phase; (v) The orbital 
analog of a breached-pair phase; (vi) The yrast effect sim¬ 
ilar to that occurring in nuclei; (vii) The temperature vs 
frequency phase diagram in a finite-size system. 

Previous work on vortices in Fermi gases considered a 
single vortex with cylindrical symmetry at zero temper¬ 
ature, either within a BdG approach in the BCS regime 
m and across the BCS-BEC crossover |551|^, or in¬ 
cluding energy-density-functional corrections at unitarity 
[211. Vortex lattices were considered only for a strictly 
two-dimensional geometry, at zero temperature in the 
BCS regime within a static [25] and a dynamic |2S] BdG 
approach, and close to at unitarity within a Ginzburg- 
Landau approach m- None of these works could address 
a comparison with the experimental data of refs. naiii], 
nor bring out the physical features highlighted above. 

We first consider the conditions of the experiment of 
ref.[T7] with N = 2 x 10® and an aspect ratio 2.5 (see 
Methods), with an atomic cloud large enough to support 
many vortices at low temperature. In this experiment, 
the trap was set in rotation at an optimal stirring fre¬ 
quency t'opt = 45Hz. 

Figure [l] shows false colour images of the array of vor¬ 
tices obtained theoretically at unitarity and zero tem¬ 
perature for the experimental set up of ref. m with D = 
f^opt = 0.8Dr, as seen from the top at z = 0 [panel (a)] 
and from the side at a: = 0 [panel (b)] (where = 27r^'r). 
In Fig. the array contains 137 vortices, whose loca¬ 
tions are not fixed a priori but are self-consistently deter¬ 
mined by solving the LPDA equation (see Supplementary 
Information). The spacing among vortices in the trian¬ 
gular array increases away from the trap center. Near the 
trap center this spacing is found to be [7i7r/(-\/3mD)]^/^, 
consistent with the Feynman’s theorem for the number 
of vortices per unit area = 2mD/7r?i which applies to 
an infinite array of vortices [2S] (where Ti is the Planck 
constant h divided by 27r). In Fig. [^, 11 vortex fila¬ 
ments are distinguished and appear to bend away from 
the z (axial) direction upon approaching the boundary 
of the cloud, so as to meet the surface perpendicularly 
with no current leaking out of the cloud (this feature 
was noted previously for a single vortex in an elongated 



FIG. 1. Arrays of vortices at unitarity for f2 = 0.8 fir- 

Magnitude of the gap parameter obtained from the self- 
consistent solution of the LPDA equation at T = 0, seen from 
the top (a) and from the side (b). c. Radial profile of the 
(magnitude of the) gap parameter (in units of its value Aq in 
between two adjacent vortices) inside the central vortex in the 
trap for various angular frequencies (in units of fir) at T = 0. 
The inset compares the vortex count over an increasing su¬ 
perfluid section of the cloud with the prediction of Feyman’s 
theorem (where Rf = ^Ef /is the Thomas-Fermi 
radius, with the Fermi energy Ef and fio given in Methods), 
d. Profile of the (magnitude of the) gap parameter (broken 
lines) and of the density (full lines) for different temperatures 
(in units of the critical temperature Tc^ for fi = 0). 


trapped bosonic cloud [29] )• Figure shows the radial 
profile of the (magnitude of the) gap parameter inside 
the central vortex for several values of the angular fre¬ 
quency D, whose shape appears to be unmodified once 
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normalized to the “asymptotic” value Ag in between two 
adjacent vortices and the radial distance p = y/x'^ + y'^ 
is expressed in terms of the local Fermi wave vector 
A:f( 0) = [37r^n(0)]^/^ where n(0) is the density at the 
trap center. [Specifically, Ag = (0.87,0.83,0.77,0.68) FIf 
for Q = (0.0, 0.4, 0.6,0.8) This is because the healing 
length of the central vortex (which is about 2 A:f( 0)“^ in 
Fig. &) is much smaller than the distance between two 
adjacent vortices (which is about 24 /cf( 0)“^ in Fig. 
where the vortex density is highest). The inset of Fig. 

shows the p-dependence of the ratio between the num¬ 
ber of vortices Afv (p) obtained numerically within a circle 
(at z = 0) with center at a; = p = 0 and radius p, and 
the corresponding number of vortices Mf{p) = n^irp^ 
expected from Feynman’s theorem, for the frequencies 
n = (0.4,0.6,0.8) Sir- Each plot terminates at the bound¬ 
ary Rs of the superfluid portion of the cloud (cf. Fig. 

which is always smaller than Rp- One concludes 
that Feynman’s theorem is satisfied, in practice, up to 
about Rs/‘2. Finally, Fig. shows the profiles of the 
(magnitude of the) gap parameter and density a,t z = 0 
along X for O = 0.8 and various temperatures in the 
superfluid phase. A non-negligible portion of the outer 
normal component of the cloud past Rs becomes normal 
due to rotation, with the result that a bi-modal distri¬ 
bution for the density emerges below a certain tempera¬ 
ture. This feature bears analogies with what occurs for 
a bosonic system, for which the superfluid component 
emerges near the trap center as a narrow peak out of a 
broad thermal distribution [30] . Our results suggest that 
the emergence of a superfluid component below could 
be detected also for a rotating fermion system directly 
from the in situ density profiles, avoiding the need to 
expand the cloud when looking for vortex arrays. 

When approaching the BCS limit or for increasing T, 
the outer normal component increases in size at the ex¬ 
pense of the superfluid component in the inner part of 
the cloud. In particular, at T = 0 the mechanism for 
activating the outer normal component stems from the 
presence of the “classical” velocity field Vn(r) = x r in 
the argument of the Fermi functions entering the coeffi¬ 
cients of the LPDA equation (see Methods). This outer 
normal component can accordingly be referred to as an 
orbital breached pair (OBP) phase, in analogy with the 
breached pair phase for imbalanced spin populations m- 
At a given temperature, the spatial extent of the OBP 
phase increases upon approaching the BCS limit as the 
magnitude |A(r)| of the gap parameter becomes progres¬ 
sively smaller, in such a way that the total number of vor¬ 
tices drops considerably when approaching this limit. [A 
related study of pair-breaking effects in a rotating Fermi 
gas along the BCS-BEC crossover was made in ref. |32j . 
albeit in the absence of vortices.] 

We have performed our calculations over an extended 
range of {kpaF)~^ across unitarity for several values 
of the angular frequency fl, and obtained the num- 



FIG. 2. Comparing the number of vortices with the ex¬ 
perimental data. a. Number of vortices Afv obtained from 
the solution of the LPDA equation over an extended coupling 
range across unitarity, for different values of the angular fre¬ 
quency D (in units of Dr) and of the temperature (in units 
of the Fermi temperature Tp = Ep/kB, fcs being the Boltz¬ 
mann constant), b. The experimental values for AC from 
ref.|17| multiplied by a factor of four (stars) are compared 
with the results of the present calculation for D = 0.8 Dr and 
two different temperatures (in units of Tp). 


ber of vortices AC shown in Fig. it at T = 0 and 
T = O.lTp. We have obtained AC from the total circula¬ 
tion ^ d£ ■ W(p = 27rAC, where the line integral is over a 
circle which encompasses the whole superfluid portion of 
the cloud at z = 0. For increasing D, the maximum of AC 
is found to shift toward the BEC side of unitarity. Past 
this maximum, the temperature has only a minor effect 
on AC- In particular at zero temperature, we find the 
maximum to occur at {kpap)~^ = (0.13,0.27,0.47) for 
D = (0.4,0.6,0.8) Dr. In all cases, a rapid decrease of AC 
occurs when approaching the BCS side of unitarity, in ac¬ 
cordance with the above argument for the presence of an 
OBP phase in the outer portion of the cloud. Figure]^ 
compares our results with the experimental values of AC 
from ref. m at the optimal stirring frequency D = 0.8 D^. 
A remarkable agreement is found, provided one multiplies 
all the experimental data by the same factor of four, ir¬ 
respective of coupling. [We attribute the discrepancy on 
the number of vortices at coupling 3.8 to the experimen¬ 
tal problems that arise deep in the BEC regime with the 
relaxation of the highest vibrational state of the molecu¬ 
lar potential involved in the Fano-Feshbach resonance of 
®Li [iT].] The need for this rescaling by a factor of four 
can be related to our finding (cf. the inset of Fig. &) 
that Feynman’s theorem is satisfied only in about one- 
fourth of the area of the cloud, leading us to speculate 
that only those vortices residing in this central portion 
of the cloud could experimentally be detected after the 
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cloud was expanded. This problem might be overcome by 
taking in situ images of the two-dimensional vortex dis¬ 
tribution, as was recently done for a Bose gas [33] . Note 
in this context also the excellent agreement with the ex¬ 
perimental data for the position of the maximum, as it 
results from Fig. for the angular frequency fl = 0.8 Vlr 
of the experiment. 

We now go on to consider the experiment of ref. [I8] . 
which measured the moment of inertia 0 of the atomic 
cloud at unitarity as a function of temperature and 
showed its progressive quenching as the temperature was 
lowered below Tc- The quenching of 0 (with respect to 
its classical value 0ci — see Methods) entails the absence 
of quantum vortices in the cloud, since these would act 
to bring the moment of inertia back to 0ci according to 
Feynman’s theorem. To get an estimate for the value 
of the frequency at which the first vortex enters the 
cloud of radial size i?s, we use the expression 

hn,,/E f = 2{kFRs)-^ HRs /£.), (i) 

where ^ is the healing length of the vortex. This ex¬ 
pression is obtained by setting E — Lft = 0, where E is 
the energy of an isolated vortex in the laboratory frame 
and L its total angular momentum [28] . The dominant 
contribution to E originates from the angular kinetic en¬ 
ergy outside the vortex core of extent With the values 
= 0.79 (from ref. [33]) and kpRs = 72. (from the 
present calculation) at unitarity and zero temperature, 
we get rici = 0.069 rir where Mir = Ep/^Q-Q is the trap 
radial angular frequency from ref. m- This value of flci 
lies within the boundaries of the full numerical calcula¬ 
tion (see Supplementary Information), which draws the 
phase diagram for the temperature dependence of the 
lower critical frequency flci about which the first vortex 
stably appears in the trap and of the upper critical fre¬ 
quency flc 2 about which the superfluid region disappears 
from the trap, in analogy to a type-II superconductor. 

Figure shows the moment of inertia <d — L/td ob¬ 
tained from the calculation of the total angular momen¬ 
tum L (see Methods) as a function of fl at T = 0 for three 
couplings across unitarity for the geometry of ref. |18j. 
The rapid rise of 0 past the coupling-dependent thresh¬ 
old flci (which decreases from the BEG to the BCS side 
of unitarity) is due to the progressive presence of vortices 
in the trap for increasing fl. In particular, the plot re¬ 
ports also the number of vortices at unitarity for a few 
values of 17, demonstrating that not too many vortices are 
needed to stabilize 0 to its classical value. The inset of 
Fig. amplihes the behavior of 0 at unitarity in a nar¬ 
row frequency range about Qci , where a smooth increase 
is found to occur before the sharp rise at flci when the 
first vortex nucleates. While this effect is essentially sup¬ 
pressed on the BEG side, the increase of 0 for 17 < 17^ 
becomes more evident toward the BGS side due to the 
progressive presence in the outer part of the cloud of the 
OBP phase discussed above. A similar effect occurs also 
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FIG. 3. Dependence of the moment of inertia on the 
angular frequency, a. Moment of inertia 0 (in units of its 
classical value 0ci - see Methods) versus the angular frequency 
fl (in units of Qr), for three characteristic couplings at zero 
temperature. At unitarity, the number of vortices entering 
the cloud are also reported for a few angular frequencies, b. 
Top view of the atomic cloud showing the array of vortices at 
unitarity for S7/17r = 0.3. 


in nuclei (where it is referred to as the yrast effect) owing 
to the finite particle number [31 [55H37] . With the total 
particle number N = 6 x 10® from ref.[T3|, this effect is 
expected to be of the order of 1/A^ « 10“® which indeed 
corresponds to the values reported in the inset of Fig. 
[^. In addition, the linear increase of 0 before I7ci ob¬ 
tained here is in line with a general argument provided 
in ref. |36| . 

In ref. HHI, the (nominal) angular frequency was es¬ 
timated at unitarity to be 0.317,., which is one order 
of magnitude larger than the threshold 17^ for the nu- 
cleation of vortices obtained by our calculation. For 
17 = 0.3 f7r our calculation predicts, in fact, that about 
20 vortices enter the cloud at unitarity and zero tem¬ 
perature, as shown in Fig. In ref.[TB], on the other 
hand, the nucleation of vortices was not considered to oc- 
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FIG. 4. Comparison with the experimental data for 
the moment of inertia at unitarity. The moment of iner¬ 
tia 0 (in units of the classical value 0ci) vs the temperature 
T (in units of the critical temperature Tc) obtained from the 
present calculation (full line) is compared with the experi¬ 
mental data of ref.|18| (dots). Reported are also the partial 
moments of inertia contributed by the inner (dashed line) and 
outer (dotted line) regions of the cloud. The inset shows the 
temperature dependence of the radius Rs (in units of Rf) 
which sets the boundary between the inner and outer regions. 


cur on the basis of a mechanism of resonant quadrupole 
mode excitation. This absence of vortices might be con¬ 
nected with a transient configuration captured by the 
experiment, while a longer time scale would be required 
to reach the situation of full thermodynamic equilibrium 
which is assumed by our calculation where vortices nu¬ 
cleate just past rici- 

Accordingly, we compare the data of ref.[TS] to our 
calculations for the moment of inertia at unitarity for 

= 0.003 fir <C flci , in such a way that no vortex ap¬ 
pears in the trap at all temperatures below Tj,. This com¬ 
parison is reported in Fig. where 0/0ci is shown versus 
T /Tc- Even though our calculated value of (= 0.37 Tp) 
differs from that estimated in ref. [TH] (— 0.21 Tp), rescal¬ 
ing the temperature by the respective values of Tc leads 
to quite a close agreement between theory and experi¬ 
ment. This kind of rescaling is a common practice in 
condensed matter and was recently considered also for 
ultra-cold Fermi gases, when comparing theoretical pre¬ 
dictions with experimental data for the temperature de¬ 
pendence of the superfluid fraction |3S]. Reported in Fig. 

are also the partial moments of inertia contributed by 
the superfluid (inner) and normal (outer) regions of the 
atomic cloud, such that 0 = 0inner + 0outer- In the inner 
region (with |r| < Rs) the superfluid and normal com¬ 
ponents coexist with each other at a given temperature, 
while in the outer region (with |r| > Rg) only the normal 
component survives. The corresponding temperature de¬ 
pendence of Rg is shown in the inset of Fig. We have 
verified that the moment of inertia in the inner region 
where A(r) 7 ^ 0 is due only to the normal component 


which is present at non-zero temperature (cf. Eq.([^ in 
Methods). Eor 12 —0, our approach reduces to the two- 
fluid model of ref. [39] , which was used recently in ref. [40] 
to compare with the experimental data of ref.[T5j for 0 . 
In ref. |40j . however, a number of numerical approxima¬ 
tions were considered to simplify the calculation, which 
are avoided in the present approach. 

The results obtained here demonstrate how vortex ar¬ 
rays evolve from the BCS to the BEC regimes when a 
Eermi superfluid is constrained in a confined geometry, 
for increasing temperature and angular velocity or a com¬ 
bination of both. Our theoretical approach has relied on 
a novel differential equation for the spatially varying gap 
parameter, which considerably reduces the storage space 
and computational time compared to the BdG equations, 
thereby advancing the current state-of-the-art standard 
for the self-consistent generation of complex arrays of 
quantum vortices under quite broad physical conditions. 


Methods 

The coupling parameter for the BCS-BEC crossover. 

The BCS-BEC crossover is spanned by the coupling param¬ 
eter (fcFflF)”^, where ap is the scattering length of the two- 
fermion problem and kp the Fermi wave vector related to 
the (trap) Fermi energy Ep by kp/{2m) = Ep = flo (3A)^^^ 
with IIo = (flrflz)^^^. Here, m is the fermion mass, N the 
total particle number, and Qr and Hz the radial and axial 
angular frequencies, respectively. [Throughout the Methods, 
we set h equal unity.] The coupling parameter ranges from 
{kFaF)~^ < —1 characteristic of the weak-coupling (BCS) 
regime when ap < 0, to {kp aF)~^ > -|-1 characteristic of the 
strong-coupling (BEC) regime when ap > 0, across the value 
(kp aF)~^ = 0 at unitarity when Iaf] diverges. 

In ref. [ 17 ], the aspect ratio between the radial {i/r = 57 Hz) 
and axial (vz = 23 Hz) trap frequencies was chosen to be 
about 2.5, in order to maximize the number of produced vor¬ 
tices at the nominal rotation frequency f = 45Hz ~ 0.8 
In ref. [18], a larger aspect ratio of about 28 between the 
radial {ur = 680 Hz) and axial (fz = 24 Hz) trap frequen¬ 
cies was instead adopted, in order to minimize the num¬ 
ber of produced vortices at the (nominal) rotation frequency 
V = 200 Hz ~ 0.3 Vr- 

LPDA equation, density, and current. The LPDA 
equation for the gap parameter A(r) was introduced in 
ref.|20|. For the present problem of a rotating trap with an¬ 
gular velocity Cl it reads: 


-^^^A(r) =Xo(r)A(r)-b2i(r)^A(r)-iXi(r)vn(r)-VA(r) 

( 2 ) 

• dk f l-2/F(E+(kir)) m\ 

(27r)3 1 2E(k|r) j 


where 


Io(r) = 


and 


Ti{r) = - 


1 f dk ( C(k|r) 


2 J (27r)3 1 2E(k|r)3 


[l-2/F(E+(klr))] 


dfF{E+{k\r)) 

dE+(k\r) 


e(k|r) k-v,(r) 1 


i7(k|r)3 mVn(r)3 E(kjr) 


• ( 4 ) 


















6 


In these expressions, C(k|r) = ^ M where fi is the 

chemical potential and V (r) the external (trapping) potential, 
£;(k]r) = A/C(klr)2 + |A(r)|2, £+(k!r) = £;(klr) - k • v„(r) 
where Vn(r) = $1 x r is the velocity field of the normal com¬ 
ponent, and fpis) = -|- 1)”'^ the Fermi function. 

Correspondingly, the expression for the current density 
within the LPDA approach in the non-rotating frame is: 

/B(£^+(klr)) . (5) 

Here, Va(r) = V:p(r)/(2m) is the velocity field of the super- 
fluid component and n(r) the number density within LPDA: 


n(r) = 


dk 

(27r)= 




In the above expessions, ^''(k|r) = C(k|r) + |rn,Vg — mvg • Vn, 
£;''(k|r) = ^e''(k|r)2 4- |A(r)|2, and (k|r) = F;''(k|r)+ k- 
(vs(r) - v„(r)). 

Note that the expression for the current is consistent 
with that of a two-fluid model (cf., e.g., ref. [41]). This can 
be seen by expanding the argument of the Fermi function in 
Eq.([^ for small values of |va(r) — Vn(r)|, thereby identiying 
the normal density through the expression 


n„(r) = - 2 


dk k^ dfE{E(k\r)) 
(27r)3 3m clif(k|r) 


(7) 


which depends on temperature through the Fermi function. 

The total angular momentum of the system in the rotating 
trap is obtained in terms of the current density |5f via the 
expression: 


L = m / dr r X j(r). 


( 8 ) 


In the text, we have indicated by L the magnitude of L. 

In addition, the moment of inertia 0 of the system is ob¬ 
tained from L = 0 D and calculated from the expression Q 
for any value of H (until the trapping potential can no longer 
sustain the atomic cloud when set into rotation). In this re¬ 
spect, the present approach is not limited to small values of 
where the formalism of linear response can be employed 
|42 |. Under these circumstances, the moment of inertia 0 
can be compared with its “classical” counterpart 0ci defined 
as follows: 


0ci = m / drn(r) (SI x r)^ 


(9) 


where n(r) is the total number density given by Eq. . In 
particular, in the limit of vanishing angular velocity 0ci co¬ 
incides with the rigid-body value 0rig as calculated from the 
density distribution, under the assumption that the whole 
cloud (including the superfluid part) performs an extremely 
slow rigid rotation whereby Vg = Vn —>■ 0 in Eq. 
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Supplementary Information 


Numerical solution of the LPDA equation 

We describe the numerical procedure that we have 
adopted in the main text to solve the LPDA equation for 
the gap parameter A(r) (cf. Eq. (2) of the main text), 
together with the number equation N = f drn{r) (with 
n{r) given by Eq. (6) of the main text) to determine the 
chemical potential /r for hxed particle number N. 

The LPDA equation can be rewritten in the form: 

V^A(r) -I- ^^A(r) - 4mfv„(r) • VA(r) = 0, (10) 

-li(r) 

where To(r) = [m^/(7raF) -l-2io(r)] and with h = 1. The 
coefficients To(r) and Xi(r) are defined by Eqs. (3) and 
(4) of the main text. The differential equation (10) is 


solved over a finite box, with the condition that A(r) 
vanishes identically from the boundary of the box out¬ 
wards. For the trap of ref. [SI] the box width is taken 
as 2.7Rp along the x and y directions and 6.8Rf along 
the z direction; for the trap of ref. [S2] the correspond¬ 
ing widths are 1.57?^ and 22Rp. In this way, the box 
is at least 1.5 times larger than the size of the cloud in 
each direction for all rotation frequencies that we have 
considered. 


The differential equation (10) is transformed into a set 


of finite-difference equations which we schematize in vec¬ 
tor form as E(A) = 0, by discretizing it over a uniform 
spatial grid. Here, A is a vector formed by the unknown 
variables A^- = A(rj) where rj is a point on the grid, 
while the equation E) = 0 corresponds to the finite dif¬ 


ference version of Eq. (10) at position r^. For the trap of 


ref. [SI] the grid is taken 800 x 800 x 30 for the x, y, and 
z directions, respectively; while for the trap of ref. [S2] 
the grid is taken 500 x 500 x 50. One is thus left with 
solving a system of Np ~ 2 x 10^ non-linear equations 
for the Np complex variables A^, where the non-linearity 
arises from the functional dependence of the coefficients 
Io(r) and Ti(r) on A(r) (cf. Eqs. (3) and (4) of the 
main text). To solve this system we have implemented 
a quasi-Netwon method as follows. The ordinary (multi¬ 
dimensional) Newton method would imply modifying A 
as follows: 


^new ^ ^old _ J 


r-1 




( 11 ) 


where J ^ is the inverse of the Jacobian matrix J with 
matrix elements Jij = —-■ Here, the Jacobian ma¬ 

trix can be calculated quite accurately with a numerical 
effort comparable to that of evaluating E(A), because 
most of the off-diagonal elements vanish and those dif¬ 
ferent from zero are linear combinations of A.^ and A^-. 
The memory storage of the sparse matrix J is set up by 
using a compressed sparse column (CSC) format. How¬ 
ever, since the numerical inversion of J is too costly, we 







have resorted to an incomplete LU factorization [S3] and 
obtained an approximate inverse of J to be inserted in 


Eq.(ll|. It is the use of this approximate inverse of J 


that makes the procedure a quasi-Newton method in¬ 
stead of an ordinary Newton method. Specifically for 
our problem, this method proves to converge better than 
alternative versions of the quasi-Newton method (such as 
the SRI or Broyden’s methods [S4]). 

For a given trial value of we routinely perfori n 4 0 
iterations for the discretized gap A according to Eq. (111. 


We then update the chemical potential /r through a single 
step of the secant method applied to the number equa¬ 
tion N = J drn{r) at fixed A(r). With this new value 
of /i, we again repeat 40 iterations for A, and so on. In 
the presence of a large number of vortices (about one 
hundred or more), 50 steps to update /i, each followed 
by 40 iterations for A, are typically required to reach a 
satisfactory convergence. With a smaller number of vor¬ 
tices, on the other hand, these numbers can considerably 
be decreased (together with the number of points for the 
spatial grid in the xy plane). In order to speed up the cal¬ 
culation (and to make it feasible, in practice, for a large 
number of vortices), the coefficients lo(i’) and 21i(r) are 
calculated over an interpolation grid 100 x 100 x 100 in 
the variables (|A|, |A|, p, = y, — V{r)), with a logarithmic 
spacing for A. [Note that this grid is over the possible 
values of |A|, |A|, and p, and not over the physical space 
spanned by the variable r.j The values of the coefficients 
2 o(r) and 21i(r) at position r are then obtained by a tri- 
linear interpolation within a cube containing the point 
|A(r)|, |A(r)|, and p{r). In this way, the most demand¬ 
ing cases (like that shown in Fig. la of the main text) 
required 30 hours of CPU time on a standard desktop 
computer (with no parallelization of the code). 

For frequencies close to flci, where the solution with 
one vortex is almost degenerate with that without vor¬ 
tices, one needs to be particularly careful. In this case, 
we have used for the initial ansatz the product of the 
gap profile Aj^pid' ~ V{x,y,z)) within a local density 
approximation for the system in the absence of rotation 
[S5] (where V(x, y, z) is the trapping potential), with the 
function 


; R-o) ■ 




x-XQ+i{y-Yo] 


which simulates a vortex of radius centered at Rq = 
{Xq,Yo). Here, the position and radius of the vortex are 
parameters that can be varied to optimize the solution. 
In particular, the initial position of the vortex should 
be slightly displaced from the trap center, to avoid the 
system being trapped in an excited state for H < Qci ■ In 
this way, for > Hcj (but still close to HcJ the vortex 
adjusts its position and radius to the convergence values. 
For H < Hcj, on the other hand, the vortex migrates 
towards the edges of the cloud and eventually disappears. 


At larger rotation frequencies, to speed up the cal¬ 
culations we have used as initial ansatz a gap profile 
given by the above local density approximation, but 
now multiplied by a triangular lattice of vortices which 
are spaced according to Feynman’s theorem. Specifi¬ 
cally, this is implemented by multiplying the gap profile 
ATF{yi—V{x, y, z)) within local density with the function 
rif=i /(r;Rv), where /(r; Rv) is given by the expression 
( [l^ with Rv = (Av, Tv) replacing Rq = (Aq, Tq) and ^v 
replacing ^o, {(Av,Tv);v = 1, - • • ,A/’v} being the initial 
positions of A/’v vortices in a triangular lattice spaced ac¬ 
cording to Feynman’s theorem. It turns out that this 
initial configuration contains about twice the number of 
vortices of the final configuration at convergence (since 
Feynman’s theorem is progressively violated when ap¬ 
proaching the border of the cloud). Correspondingly, the 
iteration procedure evolves in such a way that a number 
of vortices progressively evaporates away from the bor¬ 
der of the cloud, until the system reaches its final equi¬ 
librium configuration. Quite generally, in the course of 
the iterations the radii of the vortex cores adjust to their 
convergence values rather quickly, while a larger number 
of iterations is required to reach convergence as far as the 
positions and number of vortices are concerned. 


Spontaneous self-assembling of vortex arrays 

The distinct power of the present method is that vor¬ 
tex arrays can be generated through the cycles of self- 
consistency at a finite rotation frequency, even when the 
initial condition for the gap profile contains (essentially) 
no vortices. In this respect, during iterations we have 
found that, quite generally, vortices enter from the edges 
of the cloud and eventually reach their equilibrium posi¬ 
tions inside the cloud. 

As a demonstration of a typical numerical simula¬ 
tion that shows the evolution through the cycles of self- 
consistency, we report here a simplified version of the cal¬ 
culations presented in the main text, where now it is the 
chemical potential /i and not the particle number N to 
be kept fixed at a given value through the cycles of self- 
consistency (in this case, we use thermodynamic value 
y = Q.7b2Ep of Fig. 3b of the main text). In addition, 
to speed up the calculation further we now utilize a spa¬ 
tial grid with 300 x 300 x 35 points instead of the denser 
one with 500 x 500 x 50 points used for the calculations 
in the main text for the trap of ref. [S2]. 

Figure SI shows a typical evolution of the gap pro¬ 
file during the cycles of self-consistency for the trap 
of ref. [S2] at unitarity, zero temperature, and = 
0.3Hj., using two different initial configurations: in the 
left column, a gap profile with no vortex but only a 
phase imprint of equilateral triangular symmetry at the 
cloud edge; in the right column, a gap profile con¬ 
taining 37 vortices as required by Feynman’s theorem. 
In order to reduce the distortion of the triangular lat¬ 
tice in the left column, the equilibrium value 0.30,. 
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FIG. 5. The evolution of the gap profile through it¬ 
erations. The evolution of the gap profile in the course of 
the iterations is shown at unitarity, T — 0, and fl = O.Sflr 
for the trap parameters of ref. [S2]. Left column (from top 
to bottom): initial configuration without vortices, and config¬ 
urations after 400, 800, and 18000 iterations. Right column 
(from top to bottom): initial configuration with 37 vortices, 
and configurations after 400, 800, and 6000 iterations. 


has been reached only asymptotically in the course of 
the iteration cycles through a suitable damped saw¬ 
tooth profile of the angular frequency (a movie show¬ 
ing the complete evolution for this case is available at 

The result is 

that these two quite different initial configurations lead 
essentially (apart from an overall rotation) to the same 
final solution at convergence with a total of 18 vortices 
(note also that the panel at the bottom of the right col¬ 


http: //bcsbec.df.unicam.it/?q=node/l 



FIG. 6. Phase diagram for the critical frequencies Qc^ 
and Qc2 vs temperatnre. Temperature dependence of the 
two critical frequencies flci and f 2 c 2 (in units of the radial trap 
frequency fir) at unitarity for the trap of ref. [S2]. The shaded 
red (blue) area corresponds to the uncertainty associated with 
flci (flc 2 )- At a given temperature, no vortex is present for 
0 < fl < flci while vortex arrays appear for fli“^ < fl < f2c2 
(shaded yellow area). The broken line extrapolates the values 
of flij^ down to the point (fl = 0,T = Teg), while the dashed 
line corresponds to Eq.(l) of the main text with appropriate 
(temperature dependent) values of ^ and Rs- The inset shows 
flci at T = 0 vs the coupling (/cfAf)”^ again obtained from 
Eq.(l) of the main text with the appropriate values of ^ and 
Rs at zero temperature. 


umn of Fig. SI coincides with Fig. 3b of the main text, 
the minor differences being ascribed to the different num¬ 
ber of points in the spatial grids). 

The n vs T phase diagram for the superfluid 
phase of a neutral trapped Fermi gas 

It is interesting to combine together the two physical 
effects which yield a finite value for the moment of in¬ 
ertia in the superfluid phase, namely, the presence of: 
(i) An array of vortices even in the absence of a normal 
component when the angular frequency increases above 
a threshold, as occurs at zero temperature (cf. Fig. 3a of 
the main text); (ii) A normal component at finite tem¬ 
perature even in the absence of vortices, as occurs for 
vanishing angular frequency (cf. Fig. 4 of the main 
text). Simultaneous consideration of both effects leads 
us to construct a phase diagram for the temperature de¬ 
pendence of the lower critical frequency flci about which 
the first vortex stably appears in the trap and of the 
upper critical frequency Hc 2 about which the superfluid 
region disappears from the trap. This phase diagram is 
the analogue of that for a homogeneous type-II super¬ 
conductor, showing the temperature dependence of the 
critical magnetic fields iJci and iJca [S6]. 

The results of this calculation are reported in Fig. S2 
at unitarity for the trap corresponding to the experiment 
of ref. [S2]. More precisely, at a given temperature we 
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have found it necessary to distinguish between a lower 
(Z) and an upper (it) value for flci and for accord¬ 
ing to the following considerations. The lower value 
corresponds to the smaller angular frequency at which 
an isolated vortex placed initially close to the trap center 
(say, at a distance i?s/10 from it) begins to be attracted 
toward the trap center in the course of the cycles of the 
self-consistent solution of the LPDA equation. The up¬ 
per value corresponds instead to the smaller value of 
n at which an isolated vortex placed initially at the edge 
Rs of the superfluid part of the cloud begins to be at¬ 
tracted toward the trap center. The ensuing uncertainty 
in the identification of flci, which is unavoidably present 
for a trap with finite size, corresponds to the shaded red 
area of Fig. S2. On the other hand, the lower value 17^2 
corresponds to the smaller value of 17 at which all vor¬ 
tices have eventually disappeared from the trap, while 
the upper value is determined by the condition that 
the gap parameter itself vanishes everywhere in the trap. 
[In this context, we have found that approaching 
the spatial width of A(r) shrinks progressively but never 
becomes smaller than the size of the ground-state wave 
function of the harmonic trap, and that from this point 
on it is the height of A(r) to decrease to zero.] The ensu¬ 
ing uncertainty in the identification of 17c2 corresponds to 
the shaded blue area of Fig. S2. For the specific trap and 
coupling conditions under which Fig. S2 was constructed, 
the values of 17ci and flc 2 with their related uncertain¬ 
ties could be identified up to the maximum temperature 
0.92 Tcg where is the critical temperature in the trap 
for 17 = 0. Finally, the shaded yellow area which extends 
from 17^“^ to 17^2 is where arrays of vortices are present 
for given values of 17 and T. 

For comparison, the dashed line in Fig. S2 corresponds 
to the approximate expression (1) of the main text, in 
which we have inserted the appropriate (temperature- 
dependent) values of for an isolated vortex taken 
from ref. [S7] and of kpRs obtained from the present 
calculation. It is remarkable that the curve obtained 
from the approximate expression (1) of the main text 
falls within the shaded red area of Fig. S2 obtained by 
the full numerical calculation. We have then used this ap¬ 
proximate expression to get an estimate for the coupling 
dependence of 17^ at zero temperature which is reported 
in the inset of Fig. S2. The values of 17^ obtained in 
this way are in line with the sharp thresholds occurring 
in Fig. 3a of the main text, where, however, the numeri¬ 
cal procedure sets these thresholds at the lower values of 
17 ^') 

Chemical potential 

Figure S3 shows the coupling dependence of the chem¬ 
ical potential ^ in the trap at zero temperature for sev¬ 
eral angular frequencies approaching the limiting value 
17 = I7r, past which the fermion cloud is no longer bound. 



FIG. 7. Effect of rotation on the chemical potential. 

Coupling dependence of the chemical potential in the trap (in 
units of Ef) at T = 0 for three different angular frequencies 
(in units of fir). The inset shows the chemical potential vs Q 
for {kpap)~^ = —1. 

The rotation affects /i more markedly on the BCS than 
on the BEC side of unitarity, the dependence becoming 
rather abrupt in the BCS limit as shown in the inset of 
Fig. S3 for {kpap)~^ = —1. 

Fluctuation corrections emerging from an 
inhomogeneons mean-field approach 

Quite generally, a mean-field calculation for an inho¬ 
mogeneous situation (of the type dealt with by the BdG 
or LPDA equations) contains contributions from what 
are referred to as fluctuation corrections in a homoge¬ 
neous situation. This is because, in an inhomogeneous 
situation, the imprint of the lowest excited states can be 
found in the ground-state wave function (as discussed, 
for instance, in ref. [S8]). As an example. Fig. S4 re¬ 
ports a comparison of the coherence (healing) length of 
the gap parameter, obtained alternatively by solving the 
BdG equations for an isolated vortex (cf. ref. [S7]) and 
by adding pairing (Gaussian) fluctuations on top of the 
homogeneous BCS mean field (cf. ref. [S9]). This com¬ 
parison shows that a BdG calculation is able to capture 
fluctuation contributions beyond mean field as far as the 
spatial variations of the gap parameter are concerned. 

In addition, as emphasized in the main text, the vortex 
profile (and thus the healing length) is not affected by the 
presence of the surrounding vortices, and the distance 
between two adjacent vortices is an order of magnitude 
larger than the healing length. Possible corrections to 
the healing length should thus have a minimal impact on 
the distribution of vortices. 

Nor even the further inclusion of pairing fluctuations 
beyond the Gaussian ones is expected to change the vor¬ 
tex profile significantly. This is shown in Fig. S5, which 
compares the vortex profiles at unitarity and zero tem¬ 
perature, obtained alternatively by the BdG/LPDA ap¬ 
proach (full line) and by the Density-Functional-Theory 
(DFT) approach of ref. [SIO] (dashed-dotted and broken 
lines) within two different parameterizations (referred to 
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T/Tc 


FIG. 8. Coherence length obtained by homogeneous 
mean field plus Gaussian fiuctuations and by inho¬ 
mogeneous mean field. Comparison of the temperature 
dependence of the coherence length ^ (in units of k^^), ob¬ 
tained alternatively by including Gaussian fluctuations on top 
of the homogeneous BCS mean field (full lines) and by a nu¬ 
merical solution of the inhomogeneous BdG equations for an 
isolated vortex embedded in an infinite superfluid (dots), for 
different values of the coupling parameter {kpap)~^. The 
temperature is in units of the superfluid critical temperature 
Tc of the homogeneous system. The results of the homoge¬ 
neous calculation have been rescaled by an overall factor of 
2/3, which takes into account the different definitions used for 
the same physical quantity by the two independent numerical 
calculations. [Figure adapted from Fig. 10 of ref. [S9].] 


as I and II). Rather remarkably, the BdG/LPDA profile 
just lies within the uncertainty of the DFT profiles in this 
important case [Sll]. 
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